import numpy as np
import pandas as pd
from numpy import linalg as LA
import math

# Visualisation libraries

## plotly
from plotly.offline import init_notebook_mode, iplot 
import plotly.graph_objs as go
import plotly.offline as py
# Graphics in retina format 
%config InlineBackend.figure_format = 'retina' 

import warnings
warnings.filterwarnings('ignore')

Numerical Solution of Ordinary Differential Equations:

• Single-step methods
  • The (Forward) Euler method-Euler-method)
  • The Backward (Implicit) Euler method-Euler-method)
  • The improved Euler method
  • The Explicit Midpoint method
  • The Implicit Midpoint method
  • Trapezoidal rule
  • Heun's method
• Second-order methods
  • Verlet integration
• Second-order methods
• Multi-step Methods
  • Runge-Kutta Methods
  • Adams methods
• Taylor’s methods

def Figs(x, y, Yexact, y_range = None):
    fig = go.Figure()
    fig.add_trace(go.Scatter(x = x, y = y, mode='lines', name='Approximation'))
    fig.add_trace(go.Scatter(x = x, y = Yexact, mode='lines', name='Exact'))
    fig['layout']['xaxis'].update(range=[x[0], x[-1]])
    if y_range != None :
        fig['layout']['yaxis'].update(range=y_range)
    fig.update_xaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
                    showgrid=True, gridwidth=1, gridcolor='Lightgray',
                    title_text = '$x$')
    fig.update_yaxes(showline=True, linewidth=1, linecolor='Lightgray', mirror=True,
                    showgrid=True, gridwidth=1, gridcolor='Lightgray',
                    title_text = '$y$')
    fig.update_xaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
    fig.update_yaxes(zeroline=True, zerolinewidth=1, zerolinecolor='Lightgray')
    fig.update_layout(legend_orientation='v', plot_bgcolor= 'white', height= 700, width= 750)
    fig.update_traces(marker_line_color= 'Black', marker_line_width=0.5, opacity=1)
    fig.show()

Single-step methods

The (Forward) Euler method

Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the (forward) Euler method is defiend as follows $$y_{n+1}=y_{n}+hf(x_{n},y_{n}).$$

def forward_Euler(f, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
    y[0]=y_init
    for i in range(N):
        y[i+1]=y[i]+h*f(x[i],y[i])
    
    return x, y

Example Conisder the initial value problem $\begin{cases}y'+2xy=xe^{-x^2},\\ y(0) = 1 \end{cases},\quad 0 \leq x \leq 1$ with exact solution $y\left(x \right) = \left(1 + \dfrac{x^{2}}{2}\right) e^{- x^{2}}$, use the forward Euler to approximate the solution of the IVP.

# f(x,y(x)):
f=lambda x, y: x*np.exp(-x**(2))-2*x*y
# f'(x,y(x)):
f1=lambda x, y: np.exp(-x**2) - 2*y - 2*(x**2)*np.exp(-x**2)
# the eact solution y(x)
y_exact= lambda x: (1+(x**2)/2)*np.exp(-x**2)
# the initial value y(0)=1
y_init=1
# the domain 
x_range=[0,1]
# the step size
h=0.01
#  Forward Euler
x,y =forward_Euler(f, y_init, x_range, h)
Yexact=y_exact(x)

# Figure
Figs(x, y, Yexact, [0.5, 1.1])

it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =forward_Euler(f, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        
data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	2.4427e-02	0.00	0.00
1	1.0e-02	2.3036e-03	10.60	1.03
2	1.0e-03	2.2904e-04	10.06	1.00
3	1.0e-04	2.2891e-05	10.01	1.00

The Backward (Implicit) Euler method

Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the backward (implicit) Euler method is defiend as follows $$y_{n+1}=y_{n}+hf(x_{n+1},y_{n+1}).$$

The improved Euler method

Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the (Forward) Euler method is defiend as follows $$y_{n+1}=y_{n}+hf(x_{n},y_{n}).$$

The improved (semi-implicit) Euler method can be also is defiend as follows \begin{align} \begin{cases} z_{n}=y_{n}+hf(x_{n},y_{n})\\ y_{n+1}=y_{n}+\dfrac{h}{2}\left(f(x_{n},y_{n})+f(x_{n+1},z_{n}) \right) \end{cases} \end{align}

def Improved_forward_Euler(f, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
    y[0]=y_init
    for i in range(N):
        z=y[i]+h*f(x[i],y[i])
        y[i+1]=y[i]+(h/2)*(f(x[i],y[i])+f(x[i+1],z))
    
    return x, y

Example: Try the previous Example, using the semi-implicit Euler method.

f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01

x,y =Improved_forward_Euler(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])

it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =Improved_forward_Euler(f, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        

data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	1.0421e-01	0.00	0.00
1	1.0e-02	1.2515e-03	83.27	1.92
2	1.0e-03	1.2745e-05	98.20	1.99
3	1.0e-04	1.2768e-07	99.82	2.00

The Explicit Midpoint method

Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ The explicit midpoint method is given by the formula $$y_{n+1}=y_{n}+hf\left(x_{n}+{\frac {h}{2}},y_{n}+{\frac {h}{2}}f(x_{n},y_{n})\right).$$

def Midpoint_method(f, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
    y[0]=y_init
    for i in range(N):
        z=y[i]+h*f(x[i],y[i])
        y[i+1]=y[i]+h*f(x[i]+h/2,y[i]+(h/2)*f(x[i],y[i]))
    
    return x, y

f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01

x,y =Midpoint_method(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])

it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =Midpoint_method(f, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        

data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	1.6278e-01	0.00	0.00
1	1.0e-02	1.8020e-03	90.34	1.96
2	1.0e-03	1.8205e-05	98.98	2.00
3	1.0e-04	1.8224e-07	99.90	2.00

The Implicit Midpoint method

Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the implicit midpoint method is given by $$y_{n+1}=y_{n}+hf\left(t_{n}+{\frac {h}{2}},{\frac {1}{2}}(y_{n}+y_{n+1})\right).$$

Trapezoidal rule

Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the trapezoidal rule is given by $$y_{{n+1}}=y_{n}+{\dfrac 12}h{\Big (}f(x_{n},y_{n})+f(x_{{n+1}},y_{{n+1}}){\Big )}.$$

Heun's method

Given the initial value problem $$y'(t)=f(x,y(x)),\qquad y(x_{0})=y_{0},$$ the Heun's method is given by $$ y_{n+1}=y_{n}+{\dfrac {h}{4}}\left(f(x_{n},y_{n})+f(x_{n}+\dfrac{2}{3}h,y_{n}+\dfrac{2}{3}f(x_n,y_n))\right). $$

def Heun_method(f, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
    y[0]=y_init
    for i in range(N):
        y[i+1]=y[i]+(h/4)*(f(x[i],y[i])+3*f(x[i]+(2*h/3),y[i]+(2*h/3)*f(x[i],y[i])))
    
    return x, y

f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01

x,y =Heun_method(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])

it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =Heun_method(f, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        

data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	1.4418e-01	0.00	0.00
1	1.0e-02	1.6195e-03	89.02	1.95
2	1.0e-03	1.6386e-05	98.83	1.99
3	1.0e-04	1.6406e-07	99.88	2.00

Multi-step Methods

The Runge–Kutta method

It is defiend as follows

$y_{n+1}=y_{n}+{\dfrac {1}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\quad \text{for }n = 0, 1, 2, 3, \ldots,$ with \begin{align} \begin{cases} k_{1}&=h\ f(t_{n},y_{n}),\\k_{2}&=h\ f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{1}}{2}}\right),\\k_{3}&=h\ f\left(t_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{2}}{2}}\right),\\k_{4}&=h\ f\left(t_{n}+h,y_{n}+k_{3}\right). \end{cases} \end{align}

def runge_kutta(f, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
    y[0]=y_init
    for i in range(N):
        k1=f(x[i],y[i]);
        k2=f(x[i]+(h/2),y[i]+(h/2)*k1);
        k3=f(x[i]+(h/2),y[i]+(h/2)*k2);
        k4=f(x[i+1],y[i]+h*k3);
        y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6
    
    return x, y

f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01

x,y =runge_kutta(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])

it=3
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =runge_kutta(f, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        

data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	1.7051e-04	0.00	0.00
1	1.0e-02	2.0347e-08	8380.19	3.92
2	1.0e-03	2.0712e-12	9823.76	3.99

Adams methods

• The two-step Adams predictor-corrector method:

\begin{align} z_{n+1}&=y_{n}+\dfrac{h}{2}\left(3f(x_{n},y_{n})-f(x_{n-1},y_{n-1})\right),\\ y_{n+1}&=y_{n}+\dfrac{h}{2} \left(f(x_{n+1},z_{n+1})+f(x_{n},y_{n})\right). \end{align}

The four-step Adams predictor-corrector method uses the four-step Adams–Bashforth and Adams-Moulton methods together:

\begin{align} z_{n+1}&=y_{n}+\dfrac{h}{24} \left(55f(x_{n},y_{n})-59f(x_{n-1},y_{n-1})+37f(x_{n-2},y_{n-2})-9f(x_{n-3},y_{n-3})\right),\\ y_{n+1}&=y_{n}+\dfrac{h}{24} \left(9f(x_{n+1},z_{n+1})+19f(x_{n},y_{n})-5f(x_{n-1},y_{n-1})+f(x_{n-2},y_{n-2})\right). \end{align}

def Adams_2step(f, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
    # the first two steps in Runge-Kutta
    x, y=runge_kutta(f, y_init, x_range, h)
        
    for i in range(1,N):
        z=y[i]+(h/2)*(3*f(x[i],y[i])-f(x[i-1],y[i-1]))
        y[i+1]=y[i]+(h/2)*(f(x[i+1],z)+f(x[i],y[i]))
    
    return x, y

f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Adams_2step(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])

it=3
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =Adams_2step(f, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        

data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	7.1255e-02	0.00	0.00
1	1.0e-02	1.1251e-03	63.33	1.80
2	1.0e-03	1.1780e-05	95.50	1.98

def Adams_4step(f, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
    # the first two steps in Runge-Kutta
    x, y=runge_kutta(f, y_init, x_range, h)
        
    for i in range(3,N):
        z=y[i]+(h/24)*(55*f(x[i],y[i])-59*f(x[i-1],y[i-1])+37*f(x[i-2],y[i-2])-9*f(x[i-3],y[i-3]))
        y[i+1]=y[i]+(h/24)*(9*f(x[i+1],z)+19*f(x[i],y[i])-5*f(x[i-1],y[i-1])+f(x[i-2],y[i-2]))
    return x, y

f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Adams_4step(f, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])

it=3
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =Adams_4step(f, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        

data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	1.4631e-04	0.00	0.00
1	1.0e-02	7.0867e-08	2064.56	3.31
2	1.0e-03	7.9936e-12	8865.40	3.95

Taylor’s methods

Taylor's Method of order $n$:

\begin{align*} y(x_{i+1}) = y(x_{i})& + h f(x_{i},y(x_{i})) + \frac{h^{2}}{2} f'(x_{i},y(x_{i})) + \ldots + \frac{h^{n}}{n!}f^{(n-1)}(x_{i},y(x_{i})) \end{align*}

def Taylor_Method(f,f1, y_init, x_range, h):
    # Inputs
    # f: the ODE y'=f(t,y(x))
    # f1: the derivative of f(t,y(x))
    # y0: the initial value
    # x_range: interval [a,b]
    # h: step size h
    
    N = int((x_range[-1] - x_range[0])/h)
    x=np.linspace(x_range[0],x_range[-1],N+1)
    y=np.zeros(N+1, dtype=float)
        
    for i in range(0,N):
        y[i+1]=y[i]+h*f(x[i],y[i])+((h**2)/2)*f1(x[i],y[i])
    return x, y

f=lambda x, y: (2*y/x) + (x**2)*np.exp(x)
y_exact= lambda x: (x**2)*(np.exp(x)-np.exp(1))
y_init=0
x_range=[1,2]
h=0.01
x,y =Taylor_Method(f, f1, y_init, x_range, h)
Yexact=y_exact(x)
# Figure
Figs(x, y, Yexact, [-5, 20])

it=4
Norm=np.zeros(it, dtype=float)
h=np.zeros(it, dtype=float)
Ratio=np.zeros(it, dtype=float)
LOG=np.zeros(it, dtype=float)

for i in range(it):
    h[i]=0.1/(10**i)
for i in range(it):
    x,y =Taylor_Method(f, f1, y_init, x_range, h[i])
    Yexact=y_exact(x)
    Norm[i]=LA.norm(y-Yexact, np.inf)
    if (i>0):
        Ratio[i]=Norm[i-1]/Norm[i]
        LOG[i]=math.log(Ratio[i],10)
        

data = pd.DataFrame({'h': h, 'Norm': Norm, 'Ratio': Ratio, 'Log': LOG})
display(data.style.format({'h': "{:.1e}", 'Norm': "{:.4e}", 'Ratio': "{:.2f}", 'Log': "{:.2f}"}))

	h	Norm	Ratio	Log
0	1.0e-01	3.7985e+00	0.00	0.00
1	1.0e-02	4.4081e-01	8.62	0.94
2	1.0e-03	4.4802e-02	9.84	0.99
3	1.0e-04	4.4876e-03	9.98	1.00

---